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Why Geophysics uses Fourier Analysis 



When earth material properties are constant in any of the cartesian variables (/ . x. y. z) then 
it is useful to Fourier transform (FT) that variable. 

In seismology, the earth does not change with time (the ocean does !) so for the earth, we 
can generally gain by Fourier transforming the time axis thereby converting time-dependent 
differential equations (hard) to algebraic equations (easier) in frequency (temporal fre- 
quency). 

In seismology, the earth generally changes rather strongly with depth, so we cannot 
usefully Fourier transform the depth z axis and we are stuck with differential equations in 
z. On the other hand, we can model a layered earth where each layer has material properties 
that are constant in Then we get analytic solutions in layers and we need to patch them 
together. 

Thirty years ago, computers were so weak that we always Fourier transformed the x 
and y coordinates. That meant that their analyses were limited to earth models in which 
velocity was horizontally layered. Today we still often Fourier transform t. x. y but not z, 
so we reduce the partial differential equations of physics to ordinary differential equations 
(ODEs). A big advantage of knowing FT theory is that it enables us to visualize physical 
behavior without us needing to use a computer. 

The Fourier transform variables are called frequencies. For each axis (t, x. y. z) we 
have a corresponding frequency (x. k x . k y . k z ). The k’s are spatial frequencies, x is the 
temporal frequency. 

The frequency is inverse to the wavelength. Question: A seismic wave from the fast 
earth goes into the slow ocean. The temporal frequency stays the same. What happens to 
the spatial frequency (inverse spatial wavelength)? 

In a layered earth, the horizonal spatial frequency is a constant function of depth. We 
will find this to be Snell’s law. 

In a spherical coordinate system or a cylindrical coordinate system, Fourier transforms 
are useless but they are closely related to “spherical harmonic functions” and Bessel trans- 
formations which play a role similar to FT. 

Our goal for these four lectures is to develop Fourier transform insights and use them 
to take observations made on the earth’s surface and “downward continue” them, to extrap- 
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olate them into the earth. This is a central tool in earth imaging. 



0.0.1 Impulse response and ODEs 



When Fourier transforms are applicable, it means the “earth response” now is the same as 
the earth response later. Switching our point of view from time to space, the applicability of 
Fourier transformation means that the “impulse response” here is the same as the impulse 
response there. An impulse is a column vector full of zeros with somewhere a one, say 
(0, 0, 1, 0, 0, ■ ■ ■)' (where the prime ()' means transpose the row into a column.) An impulse 
response is a column from the matrix 
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The impulse response is the q that comes out when the input p is an impulse. In a typical 
application, the matrix would be about 1000 x 1000 and not the simple 8x6 example 
that I show you above. Notice that each column in the matrix contains the same waveform 
( 60 , 61 , 62 ). This waveform is called the “impulse response”. The collection of impulse 
responses in Equation (0.1) defines the the convolution operation. 



Not only do the columns of the matrix contain the same impulse response, but each 
row likewise contains the same thing, and that thing is the backwards impulse response 
( 6 2 , 61 , 6 0 ). Suppose ( 6 2 , 61 , 6 0 ) were numerically equal to (1, —2, 1)/Af 2 . Then equation 
(0.1) would be like the differential equation P — 9- Equation (0.1) would be a finite- 
difference representation of a differential equation. Two important ideas are equivalent; 
either they are both true or they are both false: 



1. The columns of the matrix all hold the same impulse response. 

2. The differential equation has constant coefficients. 



The story gets more complicated when we look at the boundaries, the top and bottom few 
equations. We’ll postpone that. 



0.0.2 Z transforms 

There is another way to think about equation (0.1) which is even more basic. It does not in- 
volve physics, differential equations, or impulse responses; it merely involves polynomials. 
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(That takes me back to middle school.) Let us define three polynomials. 

P{Z) = p 0 +p 1 Z + p 2 Z 2 +p 3 Z 3 +p 4 Z 4 +p 5 Z 5 +p d Z 6 (0.2) 

B{Z) = bo + hZ + ^Z 2 (0.3) 

Q(Z) — q$ + p 4 Z P q 2 Z~ + q 3 Z 3 -\- q 4 Z 4 T q§Z 3 -\- q 3 Z 3 P q^Z 1 P q$Z 8 (0.4) 

Are you able to multiply P(Z)B(Z)1 If you are, then you can examine the coefficient 

of Z 5 . You will discover that it is exactly the fifth row of equation (0.1)! Actually it is 
the sixth row because we started from zero. For each power of Z in Q(Z) = P(Z)B(Z ) 
we get one of the rows in equation (0.1). Convolution is defined to be the operation on 
polynomial coefficients when we multiply polynomials. 



0.0.3 Frequency 

The numerical value of Z doesn’t matter. It could have any numerical value. We haven’t 
needed Z to have any particular value. It happens that real values of Z lead to what are 
called Laplace transforms and complex values of Z lead to Fourier transforms. 

Let us test some numerical values of Z. Taking Z = 1/10 we notice the earliest 
coefficient in each of the polynomials is strongly emphasized in creating the numerical 
value of the polynomial, i.e., _P(1/10) = p 0 + Pi/10 + p 2 /100 + • • •. Likewise taking 
Z = 10, the latest value is strongly emphasized. This undesirable weighting of early or 
late is avoided if we use the Fourier approach and use numerical values of Z that fulfill the 
condition \Z\ = 1. Other than Z — ±1 that forces us to use complex values of Z , but there 
are plenty of those. 

Recall the complex plane where the real axis is horizontal and the imaginary axis is 
vertical. For Fourier transforms, we are interested in complex numerical values of Z which 
have unit magnitude, namely, \Z\ = 1. Examples are Z — ±1, Z — Piox Z — (l±i)/\/2. 

The numerical value Z — 1 gives what is called the zero frequency. Evaluating P{Z = 
1) = Po + Pi + P 2 + P 3 + Pa + P 5 + P<h finds the zero-frequency component of p t . The 
value Z = — 1 gives what is called the “Nyquist frequency”. P(Z = —1) = Po — p\ Pp 2 — 
P3 T p 4 — p§ T P q. The Nyquist frequency is the highest frequency that we can represent 
with sampled time functions. If our signal p t were p t — (1, —1, +1, —1, +1, —1, +1) then 
all the terms in P(Z = —1) would add together with the same polarity so that signal has a 
strong frequency component at the Nyquist frequency. 

How about frequencies inbetween zero and Nyquist? These require us to use complex 
numbers. Consider Z — i, where i = \/—l. The signal (1, i, i 2 , i 3 , i 4 , i 5 , • • •) could be 
segregated into its real and imaginary parts. The real part is (1,0, —1,0, 1,0, ■ ■ •). Its 
wavelength is twice as long as that of the Nyquist frequency so its frequency is exactly 
half. The values for Z used by Fourier transform are Z — cos u + i sin u. 

Now we will steal parts of Jon Claerbout’s books, “Earth Soundings Analysis, Process- 
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ing versus Inversion” and “Basic Earth Imaging” which are freely available on the WWW 1 . 
To speed you along though, I trim down those chapters to their most important parts. 



1 http://sep www. Stanford, edu/sep/prof/ 




Chapter 1 

Convolution and Spectra 



Time and space are ordinarily thought of as continuous, but for the purposes of computer 
analysis we must discretize these axes. This is also called “sampling” or “digitizing.” You 
might worry that discretization is a practical evil that muddies all later theoretical analysis. 
Actually, physical concepts have representations that are exact in the world of discrete 
mathematics. 



1.1 SAMPLED DATA AND Z-TRANSFORMS 



Consider the idealized and simplified signal in Figure 1.1. To analyze such an observed 



Figure 1.1: A continuous signal 

at uniform time intervals. 
[ER] 



sampled 

cs-trivl 




signal in a computer, it is necessary to approximate it in some way by a fist of numbers. 
The usual way to do this is to evaluate or observe b(t) at a uniform spacing of points in 
time, call this discretized signal b t . For Figure 1.1, such a discrete approximation to the 
continuous function could be denoted by the vector 

b t = (...0, 0,1, 2, 0,-1, -1,0, 0,...) (U) 

Naturally, if time points were closer together, the approximation would be more accurate. 
What we have done, then, is represent a signal by an abstract n- dimensional vector. 

Another way to represent a signal is as a polynomial, where the coefficients of the 
polynomial represent the value of b t at successive times. For example, 

B(Z) = 1 + 2Z + OZ 2 — Z 3 — Z 4 (1.2) 
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This polynomial is called a “Z-transform.” What is the meaning of Z here? Z should 
not take on some numerical value; it is instead the unit-delay operator. For example, the 
coefficients of ZB(Z) — Z + 2Z 2 — Z 4 — Z 5 are plotted in Figure 1.2. Figure 1.2 shows 



Figure 1.2: The coefficients of 

ZB(Z) are the shifted version of the 
coefficients of B(Z). 

[ER] 
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the same waveform as Figure 1.1, but now the waveform has been delayed. So the signal b t 
is delayed n time units by multiplying B(Z) by Z n . The delay operator Z is important in 
analyzing waves simply because waves take a certain amount of time to move from place 
to place. 

Another value of the delay operator is that it may be used to build up more complicated 
signals from simpler ones. Suppose b t represents the acoustic pressure function or the 
seismogram observed after a distant explosion. Then b t is called the “impulse response.” If 
another explosion occurred at t = 10 time units after the first, we would expect the pressure 
function y{t ) depicted in Figure 1.3. In terms of Z-transforms, this pressure function would 
be expressed as Y(Z) — B(Z) + Z 10 B(Z). 



Figure 1.3: Response to two explo- 



sions. cs-triv3 



[ER] 




1.1.1 Linear superposition 

If the first explosion were followed by an implosion of half-strength, we would have 
B(Z) — \ Z 10 B(Z ). If pulses overlapped one another in time (as would be the case if 
B(Z) had degree greater than 10), the waveforms would simply add together in the region 
of overlap. The supposition that they would just add together without any interaction is 
called the “linearity” property. In seismology we find that — although the earth is a hetero- 
geneous conglomeration of rocks of different shapes and types — when seismic waves travel 
through the earth, they do not interfere with one another. They satisfy linear superposi- 
tion. The plague of nonlinearity arises from large amplitude disturbances. Nonlinearity is 
a dominating feature in hydrodynamics, where flow velocities are a noticeable fraction of 
the wave velocity. Nonlinearity is absent from reflection seismology except within a few 
meters from the source. Nonlinearity does not arise from geometrical complications in the 
propagation path. An example of two plane waves superposing is shown in Figure 1.4. 






1.1. SAMPLED DATA AND Z-TRANSFORMS 



3 



Figure 1.4: Crossing plane waves 
superposing viewed on the left as 
“wiggle traces” and on the right as 
“raster.” 



cs-super| [ER] 





1.1.2 Convolution with Z-transform 

Now suppose there was an explosion at t — 0, a half-strength implosion at Z — 1, and 
another, quarter-strength explosion at t — 3. This sequence of events determines a “source” 
time series, x t = (1, 0, \). The Z-transform of the source is X(Z) — 1 — \Z + \Z' S . 

The observed y t for this sequence of explosions and implosions through the seismometer 
has a Z-transform Y(Z), given by 

Y(Z) = 



The last equation shows polynomial multiplication as the underlying basis of time-invariant 
linear-system theory, namely that the output Y(Z) can be expressed as the input X(Z) 
times the impulse-response filter B(Z). When signal values are insignificant except in a 
“small” region on the time axis, the signals are called “wavelets.” 

1.1.3 Convolution equation and program 

What do we actually do in a computer when we multiply two Z-transforms together? The 
filter 2 + Z would be represented in a computer by the storage in memory of the coeffi- 
cients (2, 1). Likewise, for 1 — Z, the numbers (1,-1) would be stored. The polynomial 
multiplication program should take these inputs and produce the sequence (2, —1, —1). Let 
us see how the computation proceeds in a general case, say 

X(Z)B(Z ) = Y(Z) (1.4) 

(xo + x\Z + X 2 Z~ + • • •) (bo + b\Z + & 2 Z 2 ) = yo + UiZ + ILiZ 2 + • • • (1.5) 

Identifying coefficients of successive powers of Z, we get 



B(Z) - | B(Z) + A B(Z) 

H4) 

X(Z) B(Z) 



(1.3) 
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y i 


= x^bo + Xoh 
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Xobo + Xxbi + x 0 b 2 


(1.6) 
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In matrix form this looks like 
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(1.7) 



The following equation, called the “convolution equation,” carries the spirit of the group 
shown in (1.6): 

N b 

Vk ^ ^ Z'k—jbj ( 1 - 8 ) 

i = 0 

To be correct in detail when we associate equation (1.8) with the group (1.6), we should 
also assert that either the input x k vanishes before k 0 or A), must be adjusted so that the 
sum does not extend before xq. These end conditions are expressed more conveniently by 
defining j — k — i in equation (1.8) and eliminating k getting 



N b 

Vj+i = (l- 9 ) 

t=0 

A convolution program based on equation (1.9) including end effects on both ends, is 

convolve ( ) . 



# convolution: Y(Z) = X(Z) * B(Z) 

# 



subroutine convolve ( nb, bb, nx, xx, yy ) 
integer nb # number of coefficients in 

integer nx # number of coefficients in 

# number of coefficients in 
real bb(nb) # filter coefficients 

real xx(nx) # input trace 

real YY(1) # output trace 

integer ib, ix, iy, ny 

ny = nx + nb -1 

call null ( yy, ny) 
do ib= 1, nb 

do ix= 1, nx 

yy ( ix+ib-1) = yy ( ix+ib-1) 



filter 

input 

output will be nx+nb-1 



+ xx (ix) * bb (ib) 



return; end 
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This program is written in a language called Ratfor, a “rational” dialect of Fortran. It is 
similar to the Matlab language. You are not responsible for anything in this program, but, 
if you are interested, more details in the last chapter of PVI 1 , the book that I condensed this 
from. 



1.1.4 Negative time 

Notice that X(Z) and Y ( Z ) need not strictly be polynomials; they may contain both posi- 
tive and negative powers of Z, such as 

X(Z) — • • • + d T7 h Xq + X\Z + • • • ( 1 . 10 ) 

Y(Z) = ... + 00+01 +y 0 + y 1 Z + --- ( 1 . 11 ) 

The negative powers of Z in X(Z) and Y (Z) show that the data is defined before / 0. 

The effect of using negative powers of Z in the filter is different. Inspection of ( 1 . 8 ) shows 
that the output y k that occurs at time k is a linear combination of current and previous 
inputs; that is, (ay. i < k). If the filter B(Z) had included a term like b-i/Z, then the 
output y k at time k would be a linear combination of current and previous inputs and x k +i, 
an input that really has not arrived at time k. Such a filter is called a “nonrealizable” 
filter, because it could not operate in the real world where nothing can respond now to an 
excitation that has not yet occurred. However, nonrealizable filters are occasionally useful 
in computer simulations where all the data is prerecorded. 



1.2 FOURIER SUMS 



The world is filled with sines and cosines. The coordinates of a point on a spinning wheel 
are (x. y) — (cos (ut + 0), sin(c at -\ ©)), where u> is the angular frequency of revolution and 
( f> is the phase angle. The purest tones and the purest colors are sinusoidal. The movement 
of a pendulum is nearly sinusoidal, the approximation going to perfection in the limit of 
small amplitude motions. The sum of all the tones in any signal is its “spectrum.” 

Small amplitude signals are widespread in nature, from the vibrations of atoms to the 
sound vibrations we create and observe in the earth. Sound typically compresses air by a 
volume fraction of 10 -3 to 10 -6 . In water or solid, the compression is typically 10 -6 to 
10 -9 . A mathematical reason why sinusoids are so common in nature is that laws of nature 
are typically expressible as partial differential equations. Whenever the coefficients of the 
differentials (which are functions of material properties) are constant in time and space, the 
equations have exponential and sinusoidal solutions that correspond to waves propagating 
in all directions. 

1 http://sepwww.stanford.edu/sep/prof/ 
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1.2.1 Superposition of sinusoids 

Fourier analysis is built from the complex exponential 

e~ xwt = cos ut — i sin ut (1.12) 

A Fourier component of a time signal is a complex number, a sum of real and imaginary 
parts, say 

B = Re B + ilm B (1.13) 

which is attached to some frequency. Let j be an integer and Uj be a set of frequencies. 
A signal bit) can be manufactured by adding a collection of complex exponential signals, 
each complex exponential being scaled by a complex coefficient Bj, namely, 

bit ) = ]T B :) (1.14) 

j 

This manufactures a complex-valued signal. How do we arrange for b(t) to be real? We 
can throw away the imaginary part, which is like adding b(t) to its complex conjugate bit), 
and then dividing by two: 

Re b(t) = ^ ^ (Bj e~ iWjt + B :j e™*) (1.15) 

j 

In other words, for each positive u :j with amplitude Bj, we add a negative — c jj with ampli- 
tude Bj (likewise, for every negative Uj ...). The Bj are called the “frequency function,” or 
the “Fourier transform.” Loosely, the Bj are called the “spectrum,” though in formal math- 
ematics, the word “spectrum” is reserved for the product BjBj. The words “amplitude 

spectrum” universally mean ^ BjBj. 

In practice, the collection of frequencies is almost always evenly spaced. Let j be an 
integer u — j Au so that 

b(t) = J2 B J e ~ i(jAu>)t (1.16) 

j 

Representing a signal by a sum of sinusoids is technically known as “inverse Fourier trans- 
formation.” An example of this is shown in Figure 1.5. 



1.2.2 Sampled time and Nyquist frequency 

In the world of computers, time is generally mapped into integers too, say t — nAt. This is 
called “discretizing” or “sampling.” The highest possible frequency expressible on a mesh 
is (• • • . 1 . —1, +1, —1, +1, —1, • • •), which is the same as e mn . Setting e 1Wmaxt = e mn , we 
see that the maximum frequency is 



Aiax 



7 r 

At 



(1.17) 
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:LiOw frequency Time domain 




Sum of two frequencies 



Time domain 




Figure 1.5: Superposition of two sinusoids. 



cs-cosines 



[NR] 



Time is commonly given in either seconds or sample units, which are the same when At — 
1. In applications, frequency is usually expressed in cycles per second, which is the same 
as Hertz, abbreviated Hz. In computer work, frequency is usually specified in cycles per 
sample. In theoretical work, frequency is usually expressed in radians where the relation 
between radians and cycles is u — 2nf. We use radians because, otherwise, equations are 
filled with 27 t’s. When time is given in sample units, the maximum frequency has a name: 
it is the “Nyquist frequency,” which is 7 r radians or 1/2 cycle per sample. 



1.2.3 Fourier sum 

In the previous section we superposed uniformly spaced frequencies. Now we will super- 
pose delayed impulses. The frequency function of a delayed impulse at time delay t 0 is 
e tuJto . Adding some pulses yields the ‘‘Fourier sum”: 

B{u) = = J2 b ne iwnAt (1.18) 

n n 



The Fourier sum transforms the signal b t to the frequency function B(u). Time will often 
be denoted by t, even though its units are sample units instead of physical units. Thus we 
often see b t in equations like (1.18) instead of b n , resulting in an implied At — 1. 
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1.3 FOURIER AND Z-TRANSFORM 

The frequency function of a pulse at time t n — nXt is U wnAt — (e ,jjAl )"'. The factor e ,wAt 
occurs so often in applied work that it has a name: 

Z = e iujAt (1.19) 

With this Z , the pulse at time t n is compactly represented as Z n . The variable Z makes 
Fourier transforms look like polynomials, the subject of a literature called “Z-transforms.” 
The Z- transform is a variant form of the Fourier transform that is particularly useful for 
time-discretized (sampled) functions. 

From the definition (1.19), we have Z 2 = e' ;jJ ‘ 2Ak , Z 3 = e w3At , etc. Using these equiva- 
lencies, equation (1.18) becomes 

B(u) = B(u{Z)) = J2 b nZ n (1.20) 



1.3.1 Unit circle 

In this chapter, u is a real variable, so Z = e' wAt — cos to At + i sin u> At is a complex 
variable. It has unit magnitude because sin 2 + cos 2 1. As u ranges on the real axis, Z 
ranges on the unit circle \Z\ = 1. 



1.3.2 Differentiator 

A particularly interesting factor is (1 — Z), because the filter (1, —1) is like a time derivative. 
The time-derivative filter destroys zero frequency in the input signal. The zero frequency 
is (■ • • , 1, 1, 1, • • •) with a Z-transform (• • • + Z 2 + Z 3 + Z 4 + • • •). To see that the filter 

(1 — Z) destroys zero frequency, notice that (1 — Z)( h Z 2 + Z 3 + Z 4 H ) =0. More 

formally, consider output Y(Z ) = (1 — Z)X(Z) made from the filter (1 — Z) and any 
input X(Z ). Since (1 — Z) vanishes at Z = 1, then likewise Y (Z) must vanish at Z = 1. 
Vanishing at Z = 1 is vanishing at frequency uj — 0 because Z = exp(iuAt) from (1.19). 
Now we can recognize that multiplication of two functions of Z or of u is the equivalent 
of convolving the associated time functions. 



Multiplication in the frequency domain is convolution in the time domain. 

A popular mathematical abbreviation for the convolution operator is an asterisk: equa- 
tion (1.8), for example, could be denoted by y t = x t *b t . I do not disagree with asterisk 
notation, but I prefer the equivalent expression Y (Z) = X(Z)B(Z), which simultaneously 
exhibits the time domain and the frequency domain. 

The filter (1 — Z) is often called a “differentiator.” It is displayed in Figure 1.6. 
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filter(t) 



Amp(omega) 




Figure 1.6: A discrete representation of the first-derivative operator. The filter (1, —1) is 



plotted on the left, and on the right is an amplitude response, i.e., 1 1 — Z\ versus u. cs-ddt 
[NR] 



1.3.3 Gaussian examples 



The filter (1 I Z)/2 is a running average of two adjacent time points. Applying this filter 
N times yields the filter (1 + Z) N / 2 N . The coefficients of the filter (1 + Z) N are generally 
known as Pascal’s triangle. For large N the coefficients tend to a mathematical limit 
known as a Gaussian function, exp(— a(t — 1 0 ) 2 ), where a and t 0 are constants that we 
will not determine here. We will not prove it here, but this Gaussian-shaped signal has a 
Fourier transform that also has a Gaussian shape, exp(— /5cu 2 ). The Gaussian shape is often 
called a “bell shape.” Figure 1.7 shows an example for N ~ 15. Note that, except for the 
rounded ends, the bell shape seems a good fit to a triangle function. Curiously, the filter 



filter(t) 




Amp(omega) 




Figure 1.7: A Gaussian approximated by many powers of (1 + Z). 



cs-gauss 



[NR] 



(.75 + .25 Z) N also tends to the same Gaussian but with a different / () . A mathematical 
theorem says that almost any polynomial raised to the iV-th power yields a Gaussian. 

In seismology we generally fail to observe the zero frequency. Thus the idealized 
seismic pulse cannot be a Gaussian. An analytic waveform of longstanding popularity 
in seismology is the second derivative of a Gaussian, also known as a “Ricker wavelet.” 
Starting from the Gaussian and multiplying be (1 — Z) 2 = 1 — 2 Z + Z 2 produces this old, 
favorite wavelet, shown in Figure 1.8. 
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filter(t) 




Amp(omega) 




Figure 1.8: Ricker wavelet. 



cs-ricker [NR] 



1.3.4 Inverse Z-transform 

Fourier analysis is widely used in mathematics, physics, and engineering as a Fourier 
integral transformation pair: 



r+oo . 

B( u) = / b(t) e’ wt dt (1.21) 

J — 00 

/*+° O . . 

b(t) = / B(u)e~ lU}t du (1.22) 

• 7—00 

These integrals correspond to the sums we are working with here except for some minor 
details. Books in electrical engineering redefine e twt as e~’ wt . That is like switching u to 
—u. Instead, we have chosen the sign convention of physics, which is better for wave- 
propagation studies (as explained in IEI). The infinite limits on the integrals result from 
expressing the Nyquist frequency in radians/second as n/ At. Thus, as At tends to zero, 
the Fourier sum tends to the integral. When we reach equation (??) we will see that if a 
scaling divisor of 27 r is introduced into either (1.21) or (1.22), then b(t) will equal b(t). 

The Z-transform is always easy to make, but the Fourier integral could be difficult 
to perform, which is paradoxical, because the transforms are really the same. To make 
a Z-transform, we merely attach powers of Z to successive data points. When we have 
B(Z), we can refer to it either as a time function or a frequency function. If we graph the 
polynomial coefficients, then it is a time function. It is a frequency function if we evaluate 
and graph the polynomial B(Z = e™) for various frequencies u. 



EXERCISES: 

1 Let B(Z) = 1 + Z + Z 2 + Z 3 + Z 4 . Graph the coefficients of B(Z) as a function of 
the powers of Z. Graph the coefficients of [B(Z)] 2 . 

2 As u moves from zero to positive frequencies, where is Z and which way does it rotate 

around the unit circle, clockwise or counterclockwise? 
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3 Identify locations on the unit circle of the following frequencies : ( 1 ) the zero frequency, 
(2) the Nyquist frequency, (3) negative frequencies, and (4) a frequency sampled at 10 
points per wavelength. 

4 Sketch the amplitude spectrum of Figure 1.8 from 0 to 47 T. 



1.4 CORRELATION AND SPECTRA 



The spectrum of a signal is a positive function of frequency that says how much of each 
tone is present. The Fourier transform of a spectrum yields an interesting function called 
an “autocorrelation,” which measures the similarity of a signal to itself shifted. 



1.4.1 Spectra in terms of Z-transforms 

Let us look at spectra in terms of Z-transforms. Let a spectrum be denoted S(u), where 



S(u) = \B(u)\ 2 = B(u)B(u) (1.23) 

Expressing this in terms of a three-point Z-transform, we have 

S(u) = (bo + + b 2 e- i2uj )(b 0 + + b 2 e i2 “ ) (1.24) 

S(Z) = (io + ^ + ^j(b 0 + b lZ + b 2 Z 2 ) (1.25) 

S(Z) = B(^}B(Z) (1.26) 



It is interesting to multiply out the polynomial B(l/Z ) with B(Z) in order to examine the 
coefficients of S(Z): 

S(Z) = ~ — ° g + (Mo + + b^b-i) + Cbobi + bib 2 )Z + bob 2 Z~ 

S(Z) — + so + siZ + S 2 Z 2 (1-27) 

The coefficient Sk of Z k is given by 



■^k — 'y ( bjbj+k (1.28) 

i 

Equation (1.28) is the autocorrelation formula. The autocorrelation value Sk at lag 10 
is s 10 * It is a measure of the similarity of b-, with itself shifted 10 units in time. In the 
most frequently occurring case, bj, is real; then, by inspection of (1.28), we see that the 
autocorrelation coefficients are real, and 67 ,. = .s_/ c . 
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Specializing to a real time series gives 



S(Z) 


- s 0 + Si(z + T)+s 2 (z 2 + T) 


(1.29) 


S(Z(u)) 


= So + Si(e ia; + e~ iw ) + s 2 (e i2uj + e~ i2u} ) 


(1.30) 


S M 


= s 0 + 2 si cos u + 2s 2 cos 2lo 


(1.31) 


S M 


= ^2 s k COS kiO 

1c 


(1.32) 


S(u) 


tv 

— cosine transform of 67 , 


(1.33) 



This proves a classic theorem that for real- valued signals can be simply stated as follows: 



For any real signal, the cosine transform of the autocorrelation equals the magnitude 
squared of the Fourier transform. 



1.4.2 Two ways to compute a spectrum 

There are two computationally distinct methods by which we can compute a spectrum: (1) 
compute all the s *. coefficients from (1.28) and then form the cosine sum (1.32) for each 
u; and alternately, (2) evaluate B(Z) for some value of Z on the unit circle, and multiply 
the resulting number by its complex conjugate. Repeat for many values of Z on the unit 
circle. When there are more than about twenty lags, method (2) is cheaper, because the fast 
Fourier transform discussed in chapter 2 can be used. 



1.4.3 Common signals 

Figure 1.9 shows some common signals and their autocorrelations. Figure 1.10 shows 
the cosine transforms of the autocorrelations. Cosine transform takes us from time to fre- 
quency and it also takes us from frequency to time. Thus, transform pairs in Figure 1.10 
are sometimes more comprehensible if you interchange time and frequency. The various 
signals are given names in the figures, and a description of each follows: 

cos The theoretical spectrum of a sinusoid is an impulse, but the sinusoid was truncated 
(multiplied by a rectangle function). The autocorrelation is a sinusoid under a tri- 
angle, and its spectrum is a broadened impulse (which can be shown to be a narrow 
sine- squared function). 

sine The sine function is sin(oV-) / (coot). Its autocorrelation is another sine function, and 
its spectrum is a rectangle function. Here the rectangle is corrupted slightly by 
“Gibbs sidelobes,” which result from the time truncation of the original sine. 

wide box A wide rectangle function has a wide triangle function for an autocorrelation 
and a narrow sine-squared spectrum. 
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Signal 



Autocorrelations 



Figure 1.9: Common signals and one side of their autocorrelations. | cs-autocar] [ER] 



Autocorrelations 



Spectra 



Figure 1.10: Autocorrelations and their cosine transforms, i.e., the (energy) spectra of the 
common signals, cs-spectra [ER] 
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narrow box A narrow rectangle has a wide sine- squared spectrum. 

twin Two pulses. 

2 boxes Two separated narrow boxes have the spectrum of one of them, but this spectrum is 
modulated (multiplied) by a sinusoidal function of frequency, where the modulation 
frequency measures the time separation of the narrow boxes. (An oscillation seen in 
the frequency domain is sometimes called a “quefrency.”) 

comb Fine-toothed-comb functions are like rectangle functions with a lower Nyquist fre- 
quency. Coarse-toothed-comb functions have a spectrum which is a fine-toothed 
comb. 

exponential The autocorrelation of a transient exponential function is a double-sided 
exponential function. The spectrum (energy) is a Cauchy function, 1 / (u 2 +cUq). The 
curious thing about the Cauchy function is that the amplitude spectrum diminishes 
inversely with frequency to th e first power; hence, over an infinite frequency axis, the 
function has infinite integral. The sharp edge at the onset of the transient exponential 
has much high-frequency energy. 

Gauss The autocorrelation of a Gaussian function is another Gaussian, and the spectrum 
is also a Gaussian. 

random Random numbers have an autocorrelation that is an impulse surrounded by some 
short grass. The spectrum is positive random numbers. 

smoothed random Smoothed random numbers are much the same as random numbers, 
but their spectral bandwidth is limited. 



1.4.4 Spectra of complex-valued signals 

The spectrum of a signal is the magnitude squared of the Fourier transform of the function. 
Consider the real signal that is a delayed impulse. Its Z-transform is simply Z\ so the real 
part is cos u, and the imaginary part is sin u. The real part is thus an even function of 
frequency and the imaginary part an odd function of frequency. This is also true of Z 2 and 
any sum of powers (weighted by real numbers), and thus it is true of any time function. For 
any real signal, therefore, the Fourier transform has an even real part RE and an imaginary 
odd part 10. Taking the squared magnitude gives (RE+iIO)(RE— iIO)= (RE) 2 + (IO) 2 . The 
square of an even function is obviously even, and the square of an odd function is also even. 
Thus, because the spectrum of a real-time function is even, its values at plus frequencies 
are the same as its values at minus frequencies. In other words, no special meaning should 
be attached to negative frequencies. This is not so of complex- valued signals. 

Although most signals which arise in applications are real signals, a discussion of cor- 
relation and spectra is not mathematically complete without considering complex-valued 
signals. Furthermore, complex- valued signals arise in many different contexts. In seismol- 
ogy, they arise in imaging studies when the space axis is Fourier transformed, i.e., when 
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a two-dimensional function p(t,x ) is Fourier transformed over space to P(t , fc x ). More 
generally, complex-valued signals arise where rotation occurs. For example, consider two 
vector-component wind-speed indicators: one pointing north, recording n t , and the other 
pointing west, recording w t . Now, if we make a complex- valued time series v t — n t + iw t , 
the magnitude and phase angle of the complex numbers have an obvious physical interpre- 
tation: +u> corresponds to rotation in one direction (counterclockwise), and (— uS) to rota- 
tion in the other direction. To see why, suppose n t - cos(u 0 t ; <fi) and w t = — sin^o t +</)). 
Then v t — e _ ' (aJof+0) . The Fourier transform is 

/ +OO 

e ~i(u 0 ( 1 . 34 ) 

-00 

The integrand oscillates and averages out to zero, except for the frequency u = u 0 . So the 
frequency function is a pulse at u> = u 0 : 

V{u ) = 5(uj — uj 0 )e~ l(/> (1.35) 

Conversely, if w t were sin(cj 0 ^ + <t>), then the frequency function would be a pulse at — u> 0 , 
meaning that the wind velocity vector is rotating the other way. 



1.4.5 Time-domain conjugate 

A complex-valued signal such as e lUot can be imagined as a corkscrew, where the real 
and imaginary parts are plotted on the x- and y- axes, and time t runs down the axis of the 
screw. The complex conjugate of this signal reverses the y-axis and gives the screw an 
opposite handedness. In Z-transform notation, the time-domain conjugate is written 

B(Z) = F 0 + he iu} + Ne i2u} + ■ ■ ■ (1.36) 

Now consider the complex conjugate of a frequency function. In Z-transform notation this 
is written 

5R = = bo + b[e~ iuj + l) 2 e~ i2uj + • • • (1.37) 

To see that it makes a difference in which domain we take a conjugate, contrast the two 
equations (1.36) and (1.37). The function B(^)B(Z) is a spectrum, whereas the function 
b t b t is called an “envelope function.” 

You might be tempted to think that Z — 1 /Z, but that is true only if u is real, and often 
it is not. 



1.4.6 Spectral transfer function 

Filters are often used to change the spectra of given data. With input X(Z), filters B(Z), 
and output Y(Z), we have Y(Z) = B(Z)X(Z) and the Fourier conjugate Y(l/Z) — 
B(1/Z)X(1/Z). Multiplying these two relations together, we get 

YY = (BB)(XX) (1.38) 
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which says that the spectrum of the input times the spectrum of the filter equals the spec- 
trum of the output. Filters are often characterized by the shape of their spectra; this shape 
is the same as the spectral ratio of the output over the input: 



BB 



YY 

Tx 



(1.39) 



EXERCISES: 

1 Suppose a wavelet is made up of complex numbers. Is the autocorrelation relation 
Sk = true? Is s k real or complex? Is S(ui) real or complex? 




Chapter 2 

Discrete Fourier transform 



Happily, Fourier sums are exactly invertible: given the output, the input can be quickly 
found. Because signals can be transformed to the frequency domain, manipulated there, 
and then returned to the time domain, convolution and correlation can be done faster. Time 
derivatives can also be computed with more accuracy in the frequency domain than in the 
time domain. Signals can be shifted a fraction of the time sample, and they can be shifted 
back again exactly. In this chapter we will see how many operations we associate with the 
time domain can often be done better in the frequency domain. We will also examine some 
two-dimensional Fourier transforms. 



2.1 FT AS AN INVERTIBLE MATRIX 



A Fourier sum may be written 

B{u) = ]T b t e™ = £ bt# (2.1) 

t t 

where the complex value Z is related to the real frequency lu by Z = e' w . This Fourier sum 
is a way of building a continuous function of uj from discrete signal values b t in the time 
domain. In this chapter we will study the computational tricks associated with specifying 
both time and frequency domains by a set of points. Begin with an example of a signal that 
is nonzero at four successive instants, (6 0 , &i, b 2 . b 3 ). The transform is 

B(u) = b 0 + byZ + b 2 Z 2 + b 3 Z 3 (2.2) 



The evaluation of this polynomial can be organized as a matrix times a vector, such as 



o 

cq 




'll 1 1 ' 




\ b o] 


Si 




1 w w 2 w 3 




bi 


s 2 




1 w 2 w 4 w 6 




b 2 


co 

cq 




1 w 3 w 6 w 9 




cr 

CO 



(2.3) 
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Observe that the top row of the matrix evaluates the polynomial at Z = 1, a point where 
also uj — 0. The second row evaluates B\ B(Z IT e^ 0 ), where ujq is some base 
frequency. The third row evaluates the Fourier transform for 2ujq, and the bottom row for 
3cuo- The matrix could have more than four rows for more frequencies and more columns 
for more time points. I have made the matrix square in order to show you next how we can 
find the inverse matrix. The size of the matrix in (2.3) is N = 4 . If we choose the base 
frequency Uq and hence W correctly, the inverse matrix will be 



bo 

h 


= 1/N 


'll 1 1 ' 

1 1/W 1/W 2 1/W 3 




\ B 0 ] 
B 1 




1 1/IE 2 1/IE 4 1/IE 6 




b 2 


. h _ 




1 1/IE 3 1/IE 6 1/IE 9 




. B -s _ 



Multiplying the matrix of (2.4) with that of (2.3), we first see that the diagonals are +1 as 
desired. To have the off diagonals vanish, we need various sums, such as 1 + IE + IE 2 + IE 3 
and 1 + IE 2 + if' 4 + IE 6 , to vanish. Every element (IE 6 , for example, or 1/W 9 ) is a unit 
vector in the complex plane. In order for the sums of the unit vectors to vanish, we must 
ensure that the vectors pull symmetrically away from the origin. A uniform distribution of 
directions meets this requirement. In other words, IE should be the AMh root of unity, i.e., 

W = Vl = e 2ni/N (2.5) 



The lowest frequency is zero, corresponding to the top row of (2.3). The next-to-the- 
lowest frequency we find by setting IE in (2.5) to Z = e IW °. So = 2ir/N; and for (2.4) 
to be inverse to (2.3), the frequencies required are 



Uk 



(0,1,2,..., JV — 1)27 r 
N 



(2.6) 



2.1.1 The Nyquist frequency 

The highest frequency in equation (2.6), u = 2n(N — 1 )/N, is almost 27T. This fre- 
quency is twice as high as the Nyquist frequency u = n. The Nyquist frequency is 
normally thought of as the “highest possible” frequency, because e ml , for integer t, plots 
as (• • • , 1, —1, 1, —1, 1, —1, • • •). The double Nyquist frequency function, e ?2vrt , for integer 
t, plots as (■ ■ - , 1, 1, 1, 1, 1, • • •). So this frequency above the highest frequency is really 
zero frequency! We need to recall that B(cu) = 13 (u — 2i r). Thus, all the frequencies near 
the upper end of the range (2.6) are really small negative frequencies. Negative frequen- 
cies on the interval (— 7 r, 0) were moved to interval (tt, 27t) by the matrix form of Fourier 
summation. 



2.1.2 The comb function 

Consider a constant function of time. In the frequency domain, it is an impulse at zero 
frequency. The comb function is defined to be zero at alternate time points. Multiply this 
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constant function by the comb function. The resulting signal contains equal amounts of 
two frequencies; half is zero frequency, and half is Nyquist frequency. We see this in the 
second row in Figure 2.1, where the Nyquist energy is in the middle of the frequency axis. 
In the third row, 3 out of 4 points are zeroed by another comb. We now see something like 
a new Nyquist frequency at half the Nyquist frequency visible on the second row. 



Figure 2.1: A zero-frequency 

function and its cosine transform. 
Successive rows show increasingly 
sparse sampling of the zero- 
frequency function. 

[NR] 



dft-comb 







2.1.3 Undersampled field data 

Figure 2.2 shows a recording of an airgun along with its spectrum. The original data 
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Figure 2.2: Raw data is shown on the top left, of about a half-second duration. Right shows 
amplitude spectra (magnitude of FT). In successive rows the data is sampled less densely, 
dft-undersample [ER] 



is sampled at an interval of 4 milliseconds, which is 250 times per second. Thus, the 
Nyquist frequency l/(2Al) is 125 Hz. Negative frequencies are not shown, since the 
amplitude spectrum at negative frequency is identical with that at positive frequency. Think 
of extending the top row of spectra in Figure 2.2 to range from minus 125 Hz to plus 125 
Hz. Imagine the even function of frequency centered at zero frequency — we will soon see 
it. In the second row of the plot, I decimated the data to 8 ms. This drops the Nyquist 
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frequency to 62.5 Hz. Energy that was at —10 Hz appears at 125 — 10 Hz in the second 
row spectrum. The appearance of what were formerly small negative frequencies near the 
Nyquist frequency is called “folding” of the spectrum. In the next row the data is sampled 
at 16 ms intervals, and in the last row at 32 ms intervals. The 8 ms sampling seems OK, 
whereas the 32 ms sampling looks poor. Study how the spectrum changes from one row to 
the next. 

The spectrum suffers no visible harm in the drop from 4 ms to 8 ms. The 8 ms data 
could be used to construct the original 4 ms data by transforming the 8 ms data to the 
frequency domain, replacing values at frequencies above 125/2 Hz by zero, and then inverse 
transforming to the time domain. 

(Airguns usually have a higher frequency content than we see here. Some high-frequency 
energy was removed by the recording geometry, and I also removed some when preparing 
the data.) 



2.2 INVERTIBLE SLOW FT PROGRAM 

Because Fourier sums are exactly invertible, some other things we often require can be 
done exactly by doing them in the frequency domain. 

Typically, signals are real valued. But the programs in this chapter are for complex- 
valued signals. In order to use these programs, copy the real-valued signal into a complex 
array, where the signal goes into the real part of the complex numbers; the imaginary parts 
are then automatically set to zero. 

There is no universally correct choice of scale factor in Fourier transform: choice 
of scale is a matter of convenience. Equations (2.3) and (2.4) mimic the .^-transform, 
so their scaling factors are convenient for the convolution theorem — that a product in the 
frequency domain is a convolution in the time domain. Obviously, the scaling factors of 
equations (2.3) and (2.4) will need to be interchanged for the complementary theorem that 
a convolution in the frequency domain is a product in the time domain. 



2.2.1 The slow FT code 

The s 1 owf t () routine exhibits features found in many physics and engineering programs. 
For example, the time-domain signal (which I call “tt ( ) ”), has nt values subscripted, 
from tt ( 1 ) to tt (nt ) . The first value of this signal tt ( 1 ) is located in real physical 
time at tO. The time interval between values is dt. The value of tt (it) is at time 
1 0 + (it-1 ) *dt. I do not use “if” as a pointer on the frequency axis because if is a 
keyword in most programming languages. Instead, I count along the frequency axis with a 
variable named ie. 

subroutine slowft ( adj, add, t0,dt,tt,nt, f0,df, ff,nf) 
integer it,ie, adj, add, nt, nf 
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complex cexp, cmplx, tt(nt), ff (nf ) 

real pi2, freq, time, scale, tO,dt, fO,df 

call adjnull ( adj, add, tt,2*nt, ff,2*nf) 

pi2 = 2. * 3.14159265; 
scale = 1 . /sqrt ( l.*nt) 
df = ( 1 . /dt ) / nf 

fO . 5/dt 



do ie = 1, nf { 
do it = 1, nt { 



freq= fO + df*(ie-l) 
time= tO + dt*(it-l) 



if ( adj == 0 ) 

f f (ie) = ff (ie) 

else 



tt (it) = tt (it) 



+ tt ( it ) 
+ ff (ie) 



* cexp (cmplx (0 . , pi2*freq*time) ) 

* cexp (cmplx (0 ., -pi2*freq*time) ) 



return; end 



The total frequency band is 27T radians per sample unit or 1/A t Hz. Dividing the total 
interval by the number of points nf gives A/. We could choose the frequencies to run 
from 0 to 27T radians/sample. That would work well for many applications, but it would be 
a nuisance for applications such as differentiation in the frequency domain, which require 
multiplication by — iu including the negative frequencies as well as the positive. So it 
seems more natural to begin at the most negative frequency and step forward to the most 
positive frequency. 



2.3 SYMMETRIES 



Next we examine odd/even symmetries to see how they are affected in Fourier transform. 
The even part e t of a signal b t is defined as 



The odd part is 



et = 



bt + b-t 
2 




(2.7) 



(2.8) 



By adding (2.7) and (2.8), we see that a function is the sum of its even and odd parts: 



bt — et + o t 



(2.9) 



Consider a simple, real, even signal such as (6_i, 6 0 , 6q) = (1,0,1). Its transform 
Z + 1/Z — e* w + e~ lw — 2 cos u is an even function of u, since cos u — cos(— u). 

Consider the real, odd signal (6_ x , 6 0 , bi) = (—1,0,1). Its transform Z— 1/Z = 2zsincu 
is imaginary and odd, since sinu; = — sin(— u>). 



scale 

scale 
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Likewise, the transform of the imaginary even function (i, 0, i) is the imaginary even 
function i‘2 cos u. Finally, the transform of the imaginary odd function {—i, 0, i) is real and 
odd. 

Let r and i refer to real and imaginary, e and o to even and odd, and lower-case 
and upper-case letters to time and frequency functions. A summary of the symmetries 
of Fourier transform is shown in Figure 2.3. 



Figure 2.3: Odd functions swap real 
and imaginary. Even functions do 
not get mixed up with complex num- 
bers. 



dft-reRE [NR] 




More elaborate signals can be made by adding together the three-point functions we 
have considered. Since sums of even functions are even, and so on, the diagram in Fig- 
ure 2.3 applies to all signals. An arbitrary signal is made from these four parts only, i.e., 
the function has the form b t — (re + ro) t + i(ie + io) t . On transformation of b t , each of the 
four individual parts transforms according to the table. 

Most “industry standard” methods of Fourier transform set the zero frequency as the 
first element in the vector array holding the transformed signal, as implied by equation (2.3). 
This is a little inconvenient, as we saw a few pages back. The Nyquist frequency is then the 
first point past the middle of the even-length array, and the negative frequencies lie beyond. 
Figure 2.4 shows an example of an even function as it is customarily stored. 




Figure 2.4: Even functions as customarily stored by “industry standard” FT programs. 



dft-even [NR] 
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2.3.1 Convolution in the frequency domain 

Let Y (Z) — X{Z ) B(Z). The coefficients y t can be found from the coefficients x t and b t 
by convolution in the time domain or by multiplication in the frequency domain. For the 
latter, we would evaluate both X(Z) and B(Z) at uniform locations around the unit circle, 
i.e., compute Fourier sums Xj, and Bj. from x t and b t . Then we would form Ck = XkBk 
for all k, and inverse Fourier transform to y t . The values y t come out the same as by the 
time-domain convolution method, roughly that of our calculation precision (typically four- 
byte arithmetic or about one part in 10 -6 ). The only way in which you need to be cautious 
is to use zero padding greater than the combined lengths of x t and b t . 

An example is shown in Figure 2.5. It is the result of a Fourier-domain computation 
which shows that the convolution of a rectangle function with itself gives a triangle. Notice 
that the triangle is clean — there are no unexpected end effects. 



Figure 2.5: Top shows a rectangle 
transformed to a sine. Bottom shows 
the sine squared, back transformed 
to a triangle. 

[NR] 



dft-box2triangle 






Because of the fast method of Fourier transform described next, the frequency-domain 
calculation is quicker when both X(Z) and B(Z) have more than roughly 20 coefficients. 
If either X(Z) or B(Z) has less than roughly 20 coefficients, then the time-domain calcu- 
lation is quicker. 



2.4 TWO-DIMENSIONAL FT 



Let us review some basic facts about two-dimensional Fourier transform. A two-dimensional 
function is represented in a computer as numerical values in a matrix, whereas a one- 
dimensional Fourier transform in a computer is an operation on a vector. A 2-D Fourier 
transform can be computed by a sequence of 1-D Fourier transforms. We can first trans- 
form each column vector of the matrix and then each row vector of the matrix. Alternately, 
we can first do the rows and later do the columns. This is diagrammed as follows: 

p(t, x) * — » P(t, k x ) 

I I 

P(u, x) < — > P(u, kx) 



The diagram has the notational problem that we cannot maintain the usual convention 
of using a lower-case letter for the domain of physical space and an upper-case letter for 
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the Fourier domain, because that convention cannot include the mixed objects P(t, k x ) and 
P(uj, x). Rather than invent some new notation, it seems best to let the reader rely on the 
context: the arguments of the function must help name the function. 

An example of two-dimensional Fourier transforms on typical deep-ocean data is 
shown in Figure 2.6. In the deep ocean, sediments are fine-grained and deposit slowly in 




Figure 2.6: A deep-marine dataset p(t,x) from Alaska (U.S. Geological Survey) and the 
real part of various Fourier transforms of it. Because of the long traveltime through the 
water, the time axis does not begin at t — 0. dft-plane4 [ER] 



flat, regular, horizontal beds. The lack of permeable rocks such as sandstone severely re- 
duces the potential for petroleum production from the deep ocean. The fine-grained shales 
overlay irregular, igneous, basement rocks. In the plot of P(t, k x ), the lateral continuity 
of the sediments is shown by the strong spectrum at low k x . The igneous rocks show a k x 








2.4. TWO-DIMENSIONAL FT 



25 



spectrum extending to such large k x that the deep data may be somewhat spatially aliased 
(sampled too coarsely). The plot of P(u , x) shows that the data contains no low-frequency 
energy. The dip of the sea floor shows up in (u, k x )- space as the energy crossing the origin 
at an angle. 

Altogether, the two-dimensional Fourier transform of a collection of seismograms 
involves only twice as much computation as the one-dimensional Fourier transform of each 
seismogram. This is lucky. Let us write some equations to establish that the asserted 
procedure does indeed do a 2-D Fourier transform. Say first that any function of x and t 
may be expressed as a superposition of sinusoidal functions: 

p(t,x) - j j e - iwt+ik * x P(u, kx) du dk x (2.10) 

The double integration can be nested to show that the temporal transforms are done first 
(inside): 



pit, x) 




e 7UJt P(u, k x ) du 



dk x 



j c' kx ' r P{t,k x )dk x 



The quantity in brackets is a Fourier transform over u done for each and every k x . Alter- 
nately, the nesting could be done with the k x -integral on the inside. That would imply rows 
first instead of columns (or vice versa). It is the separability of exp(— iut + i k x x) into a 
product of exponentials that makes the computation this easy and cheap. 



2.4.1 Signs in Fourier transforms 

In Fourier transforming t-, x-, and ^-coordinates, we must choose a sign convention for 
each coordinate. Of the two alternative sign conventions, electrical engineers have chosen 
one and physicists another. While both have good reasons for their choices, our circum- 
stances more closely resemble those of physicists, so we will use their convention. For the 
inverse Fourier transform, our choice is 

p(t, x , z) = j j j e - iwt + ik * x + ikzZ P(u , kx, k z ) du dk x dk z (2.1 1) 

For the forward Fourier transform, the space variables carry a negative sign, and time 
carries a positive sign. 

Let us see the reasons why electrical engineers have made the opposite choice, and why 
we go with the physicists. Essentially, engineers transform only the time axis, whereas 
physicists transform both time and space axes. Both are simplifying their lives by their 
choice of sign convention, but physicists complicate their time axis in order to simplify their 
many space axes. The engineering choice minimizes the number of minus signs associated 
with the time axis, because for engineers, d/dt is associated with iu instead of, as is the case 
for us and for physicists, with — iu. We confirm this with equation (2.11). Physicists and 
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geophysicists deal with many more independent variables than time. Besides the obvious 
three space axes are their mutual combinations, such as midpoint and offset. 

You might ask, why not make all the signs positive in equation (2.11)? The reason is 
that in that case waves would not move in a positive direction along the space axes. This 
would be especially unnatural when the space axis was a radius. Atoms, like geophysical 
sources, always radiate from a point to infinity, not the other way around. Thus, in equation 
(2.11) the sign of the spatial frequencies must be opposite that of the temporal frequency. 

The only good reason I know to choose the engineering convention is that we might 
compute with an array processor built and microcoded by engineers. Conflict of sign con- 
vention is not a problem for the programs that transform complex- valued time functions to 
complex- valued frequency functions, because there the sign convention is under the user’s 
control. But sign conflict does make a difference when we use any program that converts 
real-time functions to complex frequency functions. The way to live in both worlds is to 
imagine that the frequencies produced by such a program do not range from 0 to +7T as the 
program description says, but from 0 to — n. Alternately, we could always take the complex 
conjugate of the transform, which would swap the sign of the w-axis. 



2.4.2 Examples of 2-D FT 

An example of a two-dimensional Fourier transform of a pulse is shown in Figure 2.7. 
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Figure 2.7: A broadened pulse (left) and the real part of its FT (right). 
[ER] 



dft-ft2dofpulse 



Notice the location of the pulse. It is closer to the time axis than the frequency axis. This 
will affect the real part of the FT in a certain way (see exercises). Notice the broadening of 
the pulse. It was an impulse smoothed over time (vertically) by convolution with (1,1) and 
over space (horizontally) with (1,4,6, 4,1). This will affect the real part of the FT in another 
way. 
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Another example of a two-dimensional Fourier transform is given in Figure 2.8. This 
example simulates an impulsive air wave originating at a point on the x-axis. We see a 
wave propagating in each direction from the location of the source of the wave. In Fourier 
space there are also two lines, one for each wave. Notice that there are other lines which 
do not go through the origin; these lines are called “spatial aliases.” Each actually goes 
through the origin of another square plane that is not shown, but which we can imagine 
alongside the one shown. These other planes are periodic replicas of the one shown. 



space 
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1 /[space] 




Figure 2.8: A simulated air wave (left) and the amplitude of its FT (right). 
[ER] 



dft- airwave 



EXERCISES: 

1 Most time functions are real. Their imaginary part is zero. Show that this means that 
F(u, k ) can be determined from F{—u , —k). 

2 What would change in Figure 2.7 if the pulse were moved (a) earlier on the /-axis, and 
(b) further on the x-axis? What would change in Figure 2.7 if instead the time axis 
were smoothed with (1,4, 6, 4,1) and the space axis with (1,1)? 

3 What would Figure 2.8 look like on an earth with half the earth velocity? 

4 Numerically (or theoretically) compute the two-dimensional spectrum of a plane wave 
[S(t — px)], where the plane wave has a randomly fluctuating amplitude: say, rand(x) 
is a random number between ±1, and the randomly modulated plane wave is [(1 + 

.2 rand(x)) S(t — px)]. 

5 Explain the horizontal “layering” in Figure 2.6 in the plot of P(u, x). What determines 
the “layer” separation? What determines the “layer” slope? 
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Chapter 3 

Downward continuation of waves 



3.1 DIPPING WAVES 



We begin with equations to describe a dipping plane wave in a medium of constant veloc- 
ity. Figure 3.1 shows a ray moving down into the earth at an angle 9 from the vertical. 
Perpendicular to the ray is a wavefront. By elementary geometry the angle between the 



Figure 3.1: 
wavefront. 



Downgoin g ray 
dwnc-front [NR] 




wavefront and the earth’s surface is also 9. The ray increases its length at a speed v. The 
speed that is observable on the earth’s surface is the intercept of the wavefront with the 
earth’s surface. This speed, namely vj sin 0, is faster than v. Likewise, the speed of the 
intercept of the wavefront and the vertical axis is vj cos 9. A mathematical expression for 
a straight line like that shown to be the wavefront in Figure 3.1 is 

Z — zq — x tan 9 (3.1) 



In this expression z 0 is the intercept between the wavefront and the vertical axis. To 
make the intercept move downward, replace it by the appropriate velocity times time: 

z = —— — x tan 9 (3.2) 

cos 9 
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Solving for time gives 

t(x,z ) = - cos 8 + — sin 8 (3.3) 

v v 

Equation (3.3) tells the time that the wavefront will pass any particular location (x, z). The 
expression for a shifted waveform of arbitrary shape is f(t — t 0 ). Using (3.3) to define the 
time shift f 0 gives an expression for a wavefield that is some waveform moving on a ray. 

/ x z \ 

moving wavefield = / f — — sin 0 — - cos 8 (3.4) 

V v v ) 



3.1.1 Snell waves 

In reflection seismic surveys the velocity contrast between shallowest and deepest reflec- 
tors ordinarily exceeds a factor of two. Thus depth variation of velocity is almost always 
included in the analysis of field data. Seismological theory needs to consider waves that 
are just like plane waves except that they bend to accommodate the velocity stratification 
v{z). Figure 3.2 shows this in an idealized geometry: waves radiated from the horizontal 
flight of a supersonic airplane. The airplane passes location x at time t 0 (x) flying hori- 





Figure 3.2: Fast airplane radiating a sound wave into the earth. From the figure you can 
deduce that the horizontal speed of the wavefr ont is the same at depth z\ as it is at depth z 2 . 
This leads (in isotropic media) to Snell’s law. 



dwnc-airplane [NR] 



zontally at a constant speed. Imagine an earth of horizontal plane layers. In this model 
there is nothing to distinguish any point on the x-axis from any other point on the x-axis. 
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But the seismic velocity varies from layer to layer. There may be reflections, head waves, 
shear waves, converted waves, anisotropy, and multiple reflections. Whatever the picture 
is, it moves along with the airplane. A picture of the wavefronts near the airplane moves 
along with the airplane. The top of the picture and the bottom of the picture both move 
laterally at the same speed even if the earth velocity increases with depth. If the top and 
bottom didn’t go at the same speed, the picture would become distorted, contradicting the 
presumed symmetry of translation. This horizontal speed, or rather its inverse dt^/dx, has 
several names. In practical work it is called the stepout. In theoretical work it is called 
the ray parameter. It is very important to note that dt$/dx does not change with depth, 
even though the seismic velocity does change with depth. In a constant- velocity medium, 
the angle of a wave does not change with depth. In a stratified medium, dt 0 /dx does not 
change with depth. 

Figure 3.3 illustrates the differential geometry of the wave. The diagram shows that 




Figure 3.3: Downgoing fronts and rays in stratified medium v(z). The wavefronts are 
horizontal translations of one another, dwnc-frontz [NR] 



dt {) sin 9 

dx v 

dto cos 9 

dz v 

These two equations define two (inverse) speeds. The first is a horizontal speed, measured 
along the earth’s surface, called the horizontal phase velocity. The second is a vertical 
speed, measurable in a borehole, called the vertical phase velocity. Notice that both these 
speeds exceed the velocity v of wave propagation in the medium. Projection of wavefronts 
onto coordinate axes gives speeds larger than v, whereas projection of rays onto coordinate 
axes gives speeds smaller than v. The inverse of the phase velocities is called the stepout 
or the slowness. 

Snell’s law relates the angle of a wave in one layer with the angle in another. The 
constancy of equation (3.5) in depth is really just the statement of Snell’s law. Indeed, we 
have just derived Snell’s law. All waves in seismology propagate in a velocity-stratified 
medium. So they cannot be called plane waves. But we need a name for waves that are 
near to plane waves. A Snell wave will be defined to be the generalization of a plane 



(3.5) 

(3.6) 
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wave to a stratified medium v(z). A plane wave that happens to enter a medium of depth- 
variable velocity v(z) gets changed into a Snell wave. While a plane wave has an angle of 
propagation, a Snell wave has instead a Snell parameter p = dt 0 /dx. 

It is noteworthy that Snell’s parameter p — dt 0 /dx is directly observable at the surface, 
whereas neither v nor 9 is directly observable. Since p = dto/dx is not only observable, 
but constant in depth, it is customary to use it to eliminate 6 from equations (3.5) and (3.6): 



dt 0 

dx 


= 


sin 9 

V 


= V 


(3.7) 


dto 

dz 


= 


cos 9 

V 


II 

TT i - 1 

to 

1 

to 


(3.8) 



3.1.2 Evanescent waves 

Suppose the velocity increases to infinity at infinite depth. Then equation (3.8) tells us 
that something strange happens when we reach the depth for which p 2 equals l/v(z) 2 . 
That is the depth at which the ray turns horizontal. Below this critical depth the seismic 
wavefield damps exponentially with increasing depth. Such waves are called evanescent. 
For a physical example of an evanescent wave, forget the airplane and think about a moving 
bicycle. For a bicyclist, the slowness p is so large that it dominates 1 /v(z) 2 for all earth 
materials. The bicyclist does not radiate a wave, but produces a ground deformation that 
decreases exponentially into the earth. To radiate a wave, a source must move faster than 
the material velocity. 
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Given a vertically upcoming plane wave at the earth’s surface, say u(t,x,z = 0) = 
u(t) const (x), and an assumption that the earth’s velocity is vertically stratified, i.e. v = 
v(z), we can presume that the upcoming wave down in the earth is simply time-shifted 
from what we see on the surface. (This assumes no multiple reflections.) Time shifting 
can be represented as a linear operator in the time domain by representing it as convolution 
with an impulse function. In the frequency domain, time shifting is simply multiplying by 
a complex exponential. This is expressed as 

u(t,z ) = u(t, z — 0) * S(t + z/v) (3.9) 

U(u,z) = U(jj : z 0) (~^ z/v (3.10) 



3.2.1 Continuation of a dipping plane wave. 

Next consider a plane wave dipping at some angle 9. It is natural to imagine continuing 
such a wave back along a ray. Instead, we will continue the wave straight down. This 
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requires the assumption that the plane wave is a perfect one, namely that the same wave- 
form is observed at all x. Imagine two sensors in a vertical well bore. They should record 
the same signal except for a time shift that depends on the angle of the wave. Notice that 
the arrival time difference between sensors at two different depths is greatest for verti- 
cally propagating waves, and the time difference drops to zero for horizontally propagating 
waves. So the time shift At is u~ L cos 6 A z where 6 is the angle between the wavefront and 
the earth’s surface (or the angle between the well bore and the ray). Thus an equation to 
downward continue the wave is 

U(u : 0, z + Az) = U(uj : 9, z) exp( — iu At ) (3.11) 

( A z cos 0 \ 

— iu j (3.12) 

Equation (3.12) is a downward continuation formula for any angle 6. We can generalize 
the method to media where the velocity is a function of depth. Evidently we can apply 
equation (3.12) for each layer of thickness Az, and allow the velocity vary with z. This 
is a well known approximation that handles the timing correctly but keeps the amplitude 
constant (since |e*^| = 1) when in real life, the amplitude should vary because of reflec- 
tion and transmission coefficients. Suffice it to say that in practical earth imaging, this 
approximation is almost universally satisfactory. 

In a stratified earth, it is customary to eliminate the angle 6 which is depth variable, and 
change it to the Snell’s parameter p which is constant for all depths. Thus the downward 
continuation equation for any Snell’s parameter is 

iuAz 

v(z) 

It is natural to wonder where in real life we would encounter a Snell wave that we could 
downward continue with equation (3.13). The answer is that any wave from real life can 
be regarded as a sum of waves propagating in all angles. Thus a field data set should first 
be decomposed into Snell waves of all values of p, and then equation (3.13) can be used 
to downward continue each p, and finally the components for each p could be added. This 
process akin to Fourier analysis. We now turn to Fourier analysis as a method of downward 
continuation which is the same idea but the task of decomposing data into Snell waves 
becomes the task of decomposing data into sinusoids along the a- ax is. 




U(u,p, z + Az) — U(u,p,z) exp 



3.2.2 Downward continuation with Fourier transform 

One of the main ideas in Fourier analysis is that an impulse function (a delta function) 
can be constructed by the superposition of sinusoids (or complex exponentials). In the 
study of time series this construction is used for the impulse response of a filter. In the 
study of functions of space, it is used to make a physical point source that can manufacture 
the downgoing waves that initialize the reflection seismic experiment. Fikewise observed 
upcoming waves can be Fourier transformed over t and x. 
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Recall that a plane wave carrying an arbitrary waveform is given by equation (3.4). 
Specializing the arbitrary function to be the real part of the function exp[— iu(t — f 0 )] gives 



moving cosine wave 



cos 



u 



(— sin 0 I - cos 9 — t 
\v v 



(3.14) 



Using Fourier integrals on time functions we encounter the Fourier kernel exp(— iut). To 
use Fourier integrals on the space-axis x the spatial angular frequency must be defined. 
Since we will ultimately encounter many space axes (three for shot, three for geophone, 
also the midpoint and offset), the convention will be to use a subscript on the letter k to 
denote the axis being Fourier transformed. So k x is the angular spatial frequency on the 
rr-axis and exp (ik x x) is its Fourier kernel. For each axis and Fourier kernel there is the 
question of the sign before the i. The sign convention used here is the one used in most 
physics books, namely, the one that agrees with equation (3.14). With this convention, a 
wave moves in the positive direction along the space axes. Thus the Fourier kernel for 
(Xj z. t)- space will be taken to be 

Fourier kernel = e lkxX e tkzZ e~ iwt = exp[ i(k x x + k z z — cut)] (3.15) 



Now for the whistles, bells, and trumpets. Equating (3.14) to the real part of (3.15), 
physical angles and velocity are related to Fourier components. The Fourier kernel has the 
form of a plane wave. These relations should be memorized! 



Angles and Fourier Components 


v k x 

sin 9 = 

u 


c OS e = vk ‘ 

CU 



(3.16) 



A point in (u, k x . Ax (-space is a plane wave. The one-dimensional Fourier kernel extracts 
frequencies. The multi-dimensional Fourier kernel extracts (monochromatic) plane waves. 

Equally important is what comes next. Insert the angle definitions into the familiar 
relation sin 2 9 + cos 2 9 = 1. This gives a most important relationship: 

k 2 x + k 2 z = ^ (3.17) 

The importance of (3.17) is that it enables us to make the distinction between an arbitrary 
function and a chaotic function that actually is a wave fie Id. Imagine any function u(t. x, z). 
Fourier transform it to U (uj. k x . k z ). Look in the (uj. k x . A: 2 )-volumc for any nonvanishing 
values of U. You will have a wavefield if and only if all nonvanishing U have coordinates 
that satisfy (3.17). Even better, in practice the (t , ^(-dependence at z 0 is usually known, 
but the ^-dependence is not. Then the ^-dependence is found by assuming U is a wavefield, 
so the ^-dependence is inferred from (3.17). 

Equation (3.17) also achieves fame as the “dispersion relation of the scalar wave equa- 
tion,” a topic developed more fully in IEI. 
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3.2.3 Linking Snell waves to Fourier transforms 

To link Snell waves to Fourier transforms we merge equations (3.5) and (3.6) with equa- 
tions (3.16) 



k x 


dt 0 
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(3.18) 

(3.19) 



The basic downward continuation equation for upcoming waves in Fourier space follows 
from equation (3.13) by eliminating p by using equation (3.1 8). For analysis of real seismic 
data we introduce a minus sign because equation (3.19) refers to downgoing waves and 
observed data is made from up-coming waves. 



U (c j, k x . z + A z) = U (c o, k x , z) exp ^ \Jl - l -^f- j (3.20) 

In Fourier space we delay signals by multiplying by e' wAt , analogously, equation (3.20) 
says we downward continue signals into the earth by multiplying by e' lkzAz . Multiplication 
in the Fourier domain means convolution in time which can be depicted by the engineering 
diagram in Figure 3.4. 



input 



Filter 



P ( io , k x , 0) 




output 



P{u,k x ,z) 



Figure 3.4: Downward continuation of a downgoing wavefield. dwnc-inout [NR] 



Downward continuation is a product relationship in both the cj-domain and the k x - 
domain. Thus it is a convolution in both time and x. What does the filter look like in the 
time and space domain? It turns out like a cone, that is, it is roughly an impulse function of 
x 2 +z 2 — v 2 t 2 . More precisely, it is the Huygens secondary wave source that was exemplified 
by ocean waves entering a gap through a storm barrier. Adding up the response of multiple 
gaps in the barrier would be convolution over x. 

A nuisance of using Fourier transforms in migration and modeling is that spaces be- 
come periodic. This is demonstrated in Figure 3.5. Anywhere an event exits the frame 
at a side, top, or bottom boundary, the event immediately emerges on the opposite side. 
In practice, the unwelcome effect of periodicity is generally ameliorated by padding zeros 
around the data and the model. 
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space, kilometers 



O 0.1 0.2 0.3 0.4 0.5 





space, kilometers 
O 0.1 0.2 0.3 0.4 0.5 



Figure 3.5: A reflectivity model on the left and synthetic data using a Fourier method on 
the right. 
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3.3 A matlab program for downward continuation 

The matlab program below produces a movie. It shows a wave propating from the surface 
downward. My plan is to clean it up and make it into a little lab. 

% Wave Movies for MATLAB (thanks to JianLiang Qian) 

% : downgoing spherical wave by dispersion relation 

clear 

nt=12; nx=48; nz=96; nw=2 ; nlam=4; dx=2; dz=l; v=l; 

% allocate memory 
pp=zeros (nz , nx, nt ) ; 
q =zeros (nx, 1 ) ; 
qx=zeros (nx, 1 ) ; 
tqb= zeros (nx, nt ) ; 

bshift = zeros ( 1 , nt ) ; % column vector 

% some parameters 

lambda = nz*dz/nlam; 

pi2 = 2. *3. 141592; 

dw = v*pi2/lambda; 

dt = pi2/ (nt*dw) ; xO = nx*dx/3; zO = nz*dz/2; 
kxO = -0.5*pi2/dx; dkx = pi2/(nx*dx); 

% Generate the data for the movies 

for iw = 1 : nw % superimpose nw frequencies 

w= iw*dw; wov = w/v; % frequency / velocity 

% initial conditions for a collapsing spherical wave 

qx ( 1 : nx) =exp ( 0 . -i*wov*sqrt ( zO*zO+ ( ( 1 : nx) *dx-xO ) . * ( ( 1 : nx) *dx-xO ) ) ); 

bshift = exp ( 0 . 0-i*w*dt* (1 : nt ) ); 

for iz = l:nz % extrapolation in depth 

q=f ftshift (f ft (f ftshift (qx) ) ) ; % to Spatial freq. domain 
for ikx=2 : nx 

kx = kxO+ ( ikx-1 ) *dkx; wkx= (kx/wov) * (kx/wov) ; 
if ( wkx <= 1.0) 
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q (ikx) =exp (i*wov*sqrt (1 . O-wkx) *dz) *q ( 
end 
end 

qx =fftshift ( ifft ( fftshift ( q) ) ) ; % 

tqb= qx*bshift; % 

tqb= reshape (tqb, 1 , nx, nt ) ; % 

pp (iz, : , : ) = pp (iz, : , : ) + tqb (1, : , : ) ; % 

end 
end 



ikx) ; 



back to space domain 
an update matrix 
2d array to 3d array 
to use 3d array addition 



play the movie: (Thanks to Tapan) 
h=imagesc ( real (pp ( : , : , 1 ) ) ) ; 
colormap gray; 
set (h, ' erasemode' , ' none' ) ; 
for k=l : 6*nt 

index=mod (k, nt) +1; 

set (h, ' cdata' , real (pp ( : , : , index) ) ) ; 
drawnow; 
end 
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